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ABSTRACT 

We present a state-of-the-art TV-body code which includes a detailed treatment of 
stellar and binary evolution as well as the cluster dynamics. This code is ideal for 
investigating all aspects relating to the evolution of star clusters and their stellar pop- 
ulations. It is applicable to open and globular clusters of any age. We use the TV-body 
code to model the blue straggler population of the old open cluster M67. Preliminary 
calculations with our binary population synthesis code show that binary evolution 
alone cannot explain the observed numbers or properties of the blue stragglers. On 
the other hand, our iV-body model of M67 generates the required number of blue 
stragglers and provides formation paths for all the various types found in M67. This 
demonstrates the effectiveness of the cluster environment in modifying the nature of 
the stars it contains and highlights the importance of combining dynamics with stellar 
evolution. We also perform a series of N = 10 000 simulations in order to quantify the 
rate of escape of stars from a cluster subject to the Galactic tidal field. 

Key words: methods: numerical - stars: evolution - stars: blue stragglers - binaries: 
general - globular clusters: general - open clusters and associations: M67 



^ 1 1 INTRODUCTION 

X: 

?H , The rich environment of a star cluster provides an ideal lab- 
. . . i oratory for the study of self-gravitating systems. It also pro- 
vides important tests for stellar evolution theory and the for- 
mation of exotic stars and binaries. The colour-magnitude 
diagram (CMD) is convenient for displaying the range of 
photometrically observable stellar populations within a star 
cluster. However, in the case of dense or dynamically old 
clusters, the appearance of the CMD can be significantly 
altered by dynamical encounters between the cluster stars. 
Therefore it is necessary to combine population synthesis 
with a description of the cluster dynamics. This is impor- 
tant for both gravitational and non-gravitational interac- 
tions between stars as well as the dynamical evolution of 
the cluster as a whole. Our approach in this area is to in- 
clude a consistent treatment of stellar and binary evolution 
in a state-of-the-art iV-body code so as to allow the gener- 
ation and interaction of the full range of stellar populations 
within a cluster environment (Hurley et al. 2000). The stellar 
evolution algorithm that we use includes variable metallic- 
ity which enables us to produce realistic cluster models for 
comparison with observed cluster populations of any age. As 
the first group to have this capability we are now embarking 



on a major project to investigate the dynamical evolution of 
star clusters and their populations. 

Blue stragglers are cluster main-sequence (MS) stars 
that seem to have stayed on the main-sequence for a time ex- 
ceeding that expected from standard stellar evolution theory 
for their mass: they lie above and blueward of the turn-off 
in a cluster CMD. Ahumada & Lapasset (1995) conducted 
an extensive survey of blue straggler candidates in Galactic 
open clusters. Their results are plotted in Figure [j] as the 
mean number of blue stragglers per open cluster, relative to 
the number of main-sequence stars in the two magnitudes 
below the turn-off, as a function of the cluster age. Also 
shown in Figure ^ is the point for the old open cluster M67 
(NGC 2682) which stands out above the mean for its age, 
containing 29 proposed blue stragglers with high probabil- 
ity of cluster membership. 

Milone & Latham (1992a, hereinafter ML) have under- 
taken a long-term radial velocity observational program to 
study the blue stragglers of M67. Of a total of ten well- 
observed blue stragglers they find that six are members of 
spectroscopic binaries. One of these is a short-period binary; 
F190 with a period of 4.183 d and an eccentricity of 0.205 
(Milone & Latham 1992b). The others are long-period bi- 
naries with periods from 846 to 4913 d; three have eccentric 
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orbits and two have orbits consistent with being circular 
(Latham & Milone 1996). Thus, if the ML sample is taken 
as representative of the overall M67 blue straggler popula- 
tion, then for every ten blue stragglers we expect four to be 
single stars and six to be found in binaries with about two 
of these circular and at least four eccentric. 

This raises the question of how the blue stragglers 
formed with the most obvious scenarios involving binary 
evolution. A Case A mass transfer scenario (Kippenhahn, 
Weigert & Hoffmeister 1967) involves a main-sequence star 
filling its Roche-lobe and transferring mass to its companion, 
a less massive main-sequence star, followed by coalescence 
of the two stars as the orbit shrinks owing to angular mo- 
mentum loss. The result is a more massive main-sequence 
star that is rejuvenated relative to other stars of the same 
mass and thus evolves to become a blue straggler. This is 
an efficient method of producing single blue stragglers pro- 
vided that a large population of close binaries exists in the 
cluster. Case B mass transfer involves a main-sequence star 
accreting material from a more evolved companion and thus 
could be a likely explanation for blue stragglers in short- 
period spectroscopic binaries, such as F190. Blue stragglers 
in long-period binaries could be produced by Case C mass 
transfer when the primary is an asymptotic giant branch 
(AGB) star that has lost much of its mass. Wind accretion 
in binaries that initially have fairly large periods could also 
be responsible for such systems. 

In all these cases, except perhaps wind accretion, the 
binary orbits should be circularized by tides before and dur- 
ing mass transfer. Other scenarios are needed to explain the 
binaries in eccentric orbits and this is where the effects of 
dynamical interactions in a cluster environment become im- 
portant. Physical stellar collisions during binary- binary and 
binary-single interactions can produce blue stragglers in ec- 
centric orbits as well as allowing the possibility of an existing 
blue straggler being exchanged into an eccentric binary. Ac- 
cording to Davies (1996) encounters between binaries and 
single stars become important when the binary fraction in 
the core exceeds about 5% and binary-binary encounters 
dominate if the fraction is greater than 30%. Additionally, 
the probability of encounters depends on the density of the 
cluster. It is also possible that perturbations from passing 
stars may induce an eccentricity in a previously circular or- 
bit. As discussed by Leonard (1996) it is unlikely that any 
one formation mechanism dominates and in the case of the 
diverse blue straggler population of M67 it seems probable 
that all the above scenarios play a role. 

The aim of this paper is to compare iV-body models 
with observations of M67 to investigate the incidence and 
distribution of blue stragglers (BSs) and in so doing to con- 
strain the nature of the primordial binary population. In 
Section [| we give an overview of the observational data, in 
terms of individual stellar populations and overall cluster pa- 
rameters. We describe the details of our binary population 
synthesis in Section ^ and use it to constrain the parameters 
of the various distributions involved, with a view to maxi- 
mizing the number of blue stragglers produced. The iV-body 
code is described in Section ^ and then used in Section |^ to 
quantify the rate of escape of stars from a cluster subject 
to the tidal field of our Galaxy. In Section ^ we present our 
-/V-body model of M67 paying particular attention to the 



Log 10 (age) 

Figure 1. The number of blue stragglers relative to the num- 
ber of MS stars in the two magnitudes below the turn-off as a 
function of the population age. The stars represent the open clus- 
ter data of Ahumada & Lapasset (1995) with the M67 point an 
open symbol. The dotted line represents our population synthesis 
with a 50% binary population using the parameters of PS6 (see 
Section j^). The solid line represents our M67 iV-body simulation 
(see Section pj, note that the log-scale does not clearly indicate 
the length of the simulation). 



evolution of the blue straggler population. Our discussion of 
the results is then given, followed by conclusions. 



2 OBSERVATIONAL DATA FOR M67 

We use the CCD photometric data of Montgomery, 
Marschall & Janes (1993, hereinafter MMJ) taken from the 
Open Cluster Database (OCD: Mermilliod 1996) to con- 
struct the M67 CMD shown in Figure || Proper motion 
studies by Sanders (1977) and Girard et al. (1989) distin- 
guish stars with membership probabilities of at least 80% 
from less certain members. The CMD shows a well-defined 
photometric binary sequence, from which MMJ find that at 
least 38% of the stars in the cluster are binary systems. Of 
the BSs identified by the OCD, according to the study of 
Ahumada & Lapasset (1995), eleven are obvious candidates 
from inspection of the CMD. These eleven are all in the 
ML sample which is listed in Table ^ with each star indi- 
cated by its Sanders (1977) number. The rest of the BSs 
are much closer to the MS and its turn-off, in a clump with 
(B-V) > 0.35 and V > 11.8. The two most obvious of 
these complete the ML sample of thirteen BSs, three of 
which were rotating too rapidly to allow reliable velocity 
determinations. Of particular interest among the sample is 
the super-BS F81 which is a single star and has a mass of 
~ 3Mq (Leonard 1996), more than a factor of 2 greater 
than the cluster turn-off mass, Mto — 1.3 Mg. The binary 
S1072 is classified on the OCD as containing a BS but this 
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Figure 2. CMD for M67 (NGC 2682) using photometric data 
taken from the Open Cluster Database(OCD: Mermilliod 1996). 
Circles show probable members (P m b > 80%) as indicated by 
proper motion studies and triangles show stars of less certain 
membership. Open symbols are spectroscopic binaries. Stars iden- 
tified as blue stragglers in the OCD are plotted with stars. Open 
squares represent the seven RS CVn candidates identified by Bel- 
loni, Verbunt & Mathieu (1998, BVM). Asterisks are the six X-ray 
sources examined by van den Berg, Verbunt & Mathieu (1999, 
MVB), excluding S1082 which has already been plotted as a 
blue straggler. The cross represents the triple system observed 
by Mathieu, Latham & Griffin (1990, MLG). All the special stars 
that we have highlighted have P m b > 80%. Also plotted (full line) 
is an isochrone at t = 4 160 Myr with Z = 0.02, m-M = 9.7 and 
E{B -V) = 0.015. 



is not compatible with ubvy photometry of the system (e.g. 
van den Berg, Verbunt & Mathieu 1999, hereinafter MVB) 
so we do not count it as a BS. This leaves 28 possible BSs of 
which two have membership probability less than 80% and 
almost 2/3 have a membership probability of 98% or greater. 
Half of the 26 with P m b > 80% are the BSs that have been 
extensively studied by ML, who found that roughly half of 
their sample were in binaries. 

Belloni, Verbunt & Mathieu (1998, hereinafter BVM) 
describe the results of a ROSAT study of X-ray emission 
from M67. They detect 25 X-ray sources of which three are 
optically identified with circular binaries of orbital periods 
less than 10 d suggesting that they are RS Canum Venatico- 
rum (RS CVn) stars. These three are listed in Table |l| along 
with four other X-ray sources that are each identified with a 



spectroscopic binary (SB) that has an unknown orbital solu- 
tion and may possibly be RS CVn stars. Hall (1976) defined 
RS CVn systems to have periods between one and fourteen 
days, a hotter component of spectral type F-GV-IV, and 
strong H and K calcium lines seen in emission. The H and 
K emission is generally associated with the cool star which 
is usually a sub-giant. Multiply periodic variations in the 
light curves of these systems have been linked to spots on 
the cool star suggesting enhanced magnetic activity. This 
can be explained by rapid rotation of the cool primary star 
caused by tidal interaction with the orbit which is therefore 
likely to be circular. Mass-ratios have been determined for 
a number of eclipsing RS CVn systems (Popper 1980) and 
are generally close to one but in several cases are greater 
than one, when q — M2 /Mi . Here we let Mi denote the pri- 
mary mass and M2 the secondary mass, with the primary 
defined as the more massive star at formation of the system. 
However, the primary star is not close to filling its Roche- 
lobe in any of these systems so the mass inversion is likely 
due to a slow mass exchange such as wind accretion by the 
secondary from the primary, or simply rapid mass loss from 
the primary (Tout & Eggleton 1988). X-ray sources in M67 
which have not been identified with an optical counterpart 
by BVM are unlikely to be RS CVn systems because the 
presence of at least one sub-giant or giant would make them 
highly visible. Therefore the seven possible RS CVn systems 
is an upper limit to the expected number in M67. 

MVB obtained optical spectra for the seven X-ray 
sources found by BVM for which the X-ray emission is un- 
explained. The parameters of these systems are also listed in 
Table |l| One of these sources is the BS S1082 which was de- 
termined by ML to be single. However MVB found a second 
component in the spectrum of this star which they interpret 
as a hot sub-luminous companion. Also among the sample 
are two so-called subsubgiants whose nature is not yet un- 
derstood. 

The metallicity given for M67 on the OCD is solar al- 
though values determined by other authors would indicate 
that it is slightly sub-solar, e.g. [Fe/H] = -0.04±0.12 (Hobbs 
& Thorburn 1991) and [Fe/H] = -0.09±0.07 (Friel & Janes 
1993). The reddening of M67 is reported by various authors 
to be either E(B - V) = 0.032 (Nissen, Twarog & Craw- 
ford 1987), E(B - V) = 0.034 ± 0.019 (Fan et al. 1996) 
or E(B - V) = 0.05 (MM J), which is rather small con- 
sidering its distance of about 800 pc (OCD) from the Sun. 
M67 is an old open cluster with an age of 4 to 5 Gyr. Car- 
raro et al. (1994) used Z ~ 0.016, a distance modulus of 
m - M = 9.5 and E(B - V) = 0.02 to derive an age of 
4.8 Gyr with isochrones based on stellar models that in- 
cluded some convective overshooting. Fan et al. (1996) used 
a similar metallicity and distance modulus to obtain an age 
of 4.0 Gyr from their CCD observations of M67 while the 
OCD gives an age of 5.2 Gyr using m — M = 9.75. The de- 
tailed models of Pols et al. 1998, also including convective 
overshoot, used in conjunction with Z — 0.017, m — M = 9.6 
and E(B - V) = 0.032 give an age of 4.17 Gyr. From the 
white dwarf cooling sequence observed in M67 Richer et 
al. (1998) deduce an age of 4.3 Gyr. Shown on Figure ^ 
is an isochrone at t — 4.16 Gyr computed using the stel- 
lar evolution formulae of Hurley, Pols & Tout (2000) with 
Z = 0.02, m-M = 9.7 and E{B - V) = 0.015. We con- 
vert to observed colours with bolometric corrections com- 
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ID # 


V 


B — V 


P(d) 


e 


Rcf. 


comments 


S0752 


11.32 


0.29 


1003 


0.317 


ML 


BS binary 


S0968 


11.28 


0.13 






ML 


BS single 


S0975 


11.08 


0.43 


1221 


0.088 


ML 


BS binary 


S0977 


10.03 


-0.07 






ML 


BS single (F81) 


S0997 


12.13 


0.46 


4913 


0.342 


ML 


BS binary 


olUDO 


1 n go 


U.J.1 






IVLJj 


BS, fast rotator 


S1082 


11.25 


0.41 






ML 


BS single 


S1195 


12.34 


0.42 


1154 


0.066 


ML 


BS binary 


S1263 


11.06 


0.19 






ML 


BS single 


S1267 


10.57 


0.27 


846 


0.475 


ML 


BS binary 


S1280 


12.26 


0.26 






ML 


BS, fast rotator 


S1284 


10.94 


0.22 


4.18 


0.205 


ML 


BS binary (F190) 


S1434 


10.70 


0.11 






ML 


BS, fast rotator 


S0760 


13.29 


0.57 






BVM 


SB, poss. RS CVn 


S0972 


15.37 


0.89 






BVM 


SB, poss. RS CVn 


S0999 


12.60 


0.78 


10.06 


0.00 


BVM 


RS CVn 


S1019 


14.34 


0.81 






BVM 


SB, poss. RS CVn 


S1045 


12.54 


0.59 


7.65 


0.00 


BVM 


RSCVn 


S1070 


13.90 


0.62 


2.66 


0.00 


BVM 


RSCVn 


S1077 


12.60 


0.64 






BVM 


SB, poss. RS CVn 


S1040 


11.52 


0.87 


42.83 


0.027 


MVB 


giant + WD 


S1063 


13.79 


1.05 


18.39 


0.217 


MVB 


subsubgiant 


S1072 


11.31 


0.61 


1495 


0.32 


MVB 


BS in OCD 


S1082 


11.25 


0.41 






MVB 


BS, poss. companion? 


S1113 


13.63 


0.99 


2.823 


0.031 


MVB 


subsubgiant 


S1237 


10.78 


0.94 


697.8 


0.105 


MVB 


eccentric binary 


S1242 


12.72 


0.68 


31.78 


0.664 


MVB 


eccentric binary 


S1234 


12.65 


0.57 


4.36 


0.06 


MLG 


triple 



Table 1. Selected M67 stars identified with various observed samples. For S1234 the parameters are for the inner orbit. 



puted by Kurucz (1992) from synthetic stellar spectra. This 
isochrone gives a good match between the MS hook and 
the observed gap in stars just below the end of the MS. It 
also gives a good fit to the giant branch (GB) and the red 
clump of stars on the zero-age horizontal branch, even with 
the uncertainties in the model atmospheres for cool stars 
used to compute the colour conversions. It should however 
be noted that an isochrone with Z — 0.017, m — M = 9.75 
and E(B — V) — 0.04 gives an equally good match for the 
same age and other combinations can be found by varying 
the age. 



3 BINARY POPULATION SYNTHESIS 

To investigate the data in Figure |l| we consider a series 
of binary populations evolved from various initial distribu- 
tions of orbital separations, mass-ratios and eccentricities 
in order to determine the number of blue stragglers pro- 
duced. For each run we evolve two million binaries to an 
age of 10 10 yr with the binary population synthesis code de- 
scribed by Hurley, Tout & Pols (2000). This code, which 
represents a thorough reworking of that presented by Tout 
et al. (1997), includes tidal circularization and synchroniza- 
tion, stellar rotation, angular momentum loss mechanisms, 
common-envelope evolution and supernova kicks, in addi- 
tion to Roche-lobe overflow and mass transfer by a stellar 
wind. 



Initially we vary only the separation distribution, tak- 
ing the eccentricity as uniformly distributed. We choose the 
binary mass from the initial mass function (IMF) of Kroupa, 
Tout & Gilmore (1991, KTG1), by means of the generating 
function 



Mb 

M<T 



= 0.33 



1 



(i-x) u -'° + om(i-xy 



1.04 



,(1) 



where X is uniformly distributed between appropriate lim- 
its to give 0.2 < M b /M Q < 100. Because it has not 
been corrected for the effect of binaries it is more correctly 
used for total system mass. We then choose the component 
masses according to a uniform distribution of mass-ratio 
constrained by the single star limits of O.IMq and 50. 0M©, 
i.e. 



: (0.1/ (M b -0.1), 0.02 (M b -50)) < q < 1. 



(2) 



We follow Eggleton, Fitchett & Tout (1989, hereinafter 
EFT) by taking the distribution of orbital separations ex- 
pressed as 



(-) 



= sec (kW) + tan (kW) , 



(3) 



where W 6 [—1,1], and uniformly distributed and k satisfies 



(4) 



This distribution is symmetric in log a about a peak at a n 
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Run 


Z 


Separation 




e 


q 


-/Vbs 




PS1 


0.02 


EFT30 




uniform 


KTG1, uniform 


1.1 


1.9 


PS2 


0.02 


EFT17 




uniform 


KTG1, uniform 


1.9 


2.4 


PS3 


0.02 


EFT10 




n n lfrvrm 

U.1111W1 111 


KTOI uniform 

XV X VJ ±- L11111W1111 


3.2 


2.7 


PS4 


0.02 


FLAT LOG 




uniform 


KTG1, uniform 


8.2 


4.3 


PS5 


0.02 


EFT10, max 


200 


uniform 


KTG1, uniform 


3.9 


3.3 


PS6 


0.02 


EFT10, max 


200 


thermal 


KTG1, uniform 


4.3 


3.5 


PS7 


0.02 


EFT10, max 


200 


uniform 


KTG1, min 0.8 


3.4 


3.3 


PS8 


0.02 


EFT10, max 


200 


uniform 


KTG3 


2.8 


2.9 


PS9 


0.01 


EFT10, max 200 


thermal 


KTG1, uniform 


5.3 


3.5 



Table 2. 

4.2 Gyr. 



Numbers are per 5 000 binaries, i.e. 10 000 stars, at 



and ranges from a minimum separation of (a m to a maxi- 
mum of asm/C- We choose the constants £ and /3 to be 10~ 3 
and 0.33 respectively. A choice of a m ~ 30 AU corresponds 
to the Gaussian-like period distribution for nearby solar-like 
stars found by Duquennoy & Mayor (1991) which has a peak 
period P ~ 180 yr. 

The results are shown in Table |2| where EFTa; represents 
separations chosen from eq. (|^) with a peak, a m , at xAU. 
The distribution flat in log a used in PS4 has the same limits 
as the EFT10 distribution of run PS3. For each run the 
number of BSs present at 4.2 Gyr per 5 000 binaries is shown. 
Our BSs are MS stars that have a mass at least 2% greater 
than the cluster turn-off mass at that time, defined as the 
stellar mass which is currently due to leave the MS. Also 
shown is the number of RS CVn stars, which may provide 
an additional constraint to the population. We define these 
as circular binaries with P < 20 d containing a sub-giant or 
giant primary losing mass in a wind, some of which is being 
accreted by the secondary. This ensures that the primary 
is rotating faster than it would as a single star at the same 
evolutionary stage so that the system is magnetically active. 

We model another five populations using the EFT10 
separation distribution with an upper limit of 200 AU. The 
peak at 10 AU in the distribution is still consistent with the 
Duquennoy & Mayor (1991) observations as is the maxi- 
mum of 200 AU. This also agrees with the findings of Math- 
ieu, Latham & Griffin (1990, hereinafter MLG) for 22 spec- 
troscopic binaries in M67 while the maximum of 200 AU is 
greater than the hard/soft binary limifj^] (Heggie 1975) ex- 
pected for an open cluster. Results for runs using EFT10 
with varying choices for the eccentricity and mass-ratios are 
also given in Table ^. The IMF used in run PS8 is the sin- 
gle star IMF of Kroupa, Tout & Gilmore (1993, KTG3) and 
each binary component is chosen independently. The metal- 
licity for each run is Z — 0.02 except for PS9 which examines 
the effect of using a lower value. 

Run PS4 generates the most blue stragglers. We show 
in Section B that M67 probably contained about 40 000 stars 



* Heggie (1975) defined hard binaries to be sufficiently close that 
their binding energy exceeds the mean kinetic energy of the clus- 
ter stars. For binary component masses similar to the mean stel- 
lar mass of the cluster this roughly translates to a separation less 
than dhard = 2rj 1 /A r where is the half-mass radius and N is 
the number of cluster stars. 




* -* * *. * 



10 20 30 40 50 

period (d) 

Figure 3. Eccentricity versus Period distribution showing the 
initial conditions that evolved to BSs via Case A mass transfer 
(solid stars, single stragglers) and via Case B mass transfer (open 
stars, stragglers in close circular binaries) for the run PS6. Note 
that some wide Case C mass transfer binaries were formed at 
P > 50 d. 



initially: roughly 13 500 binaries with a 50% fraction. If all 
binaries that produce a BS are retained by the cluster then 
PS4 can explain the number found in M67. However, only 
about 25% of the BSs are in binaries, all of these have cir- 
cular orbits, and practically none are found in binaries with 
periods greater than a year. So dynamical encounters during 
cluster evolution are required to explain BSs found in wide 
binaries and those found with eccentric orbits. Although the 
observations described by Abt (1983) are consistent with a 
flat distribution of log a, more recent surveys (e.g. Duquen- 
noy & Mayor 1991) favour a peaked distribution such as is 
used in all runs except PS4. The flat distribution has also 
been ruled out by EFT. Of the runs with a peaked distribu- 
tion, PS6, with a thermal eccentricity distribution (Heggie 
1975), produces the most BSs (excluding run PS9 which 
has lower metallicity). The initial systems in run PS6 that 
lead to the formation of BSs present at 4.2 Gyr are shown 
in Figure |[ Noticeably evident is the effectiveness of tidal 
circularization at bringing eccentric binaries close enough to 
make mass transfer possible. All instances of Case A mass 
transfer lead to coalescence of the two MS stars. Of the BSs 
produced by PS6 71.9% are single, 27.8% are in close bi- 
naries resulting from Case B mass transfer and only 0.3% 
are in wide binaries produced by Case C mass transfer or 
wind accretion. Figure ^ shows the mass distribution of BSs 
present at 4.2 Gyr. Binary evolution alone cannot explain 
BSs with mass greater than twice the cluster turn-off mass, 
such as F81 in M67. 

The evolution of run PS6 is shown in Figure [j] for a 
50% binary fraction and the assumption that binaries with 
mass-ratios less than 0.4 would be observed as MS stars. 
Even though the ratio of BSs to bright MS stars is decreas- 
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binary 



2.5 



M/M- 



Figure 4. The distribution of stellar mass for blue stragglers at 
4.2 Gyr in run PS6. The distribution is shown for all blue strag- 
glers (hollow) and those that are in binaries (hatched). The cluster 
turn-off mass at this time is Mxo = 1-30 Mq. 



metallicity so that binary interaction is more likely and TV B s 
in PS9 is always greater than in PS6. 

Figure ^ shows that, as the population evolves, binary 
evolution alone cannot account for the number of observed 
BSs in open clusters if a realistic separation distribution is 
used. Cluster dynamics is therefore not only important for 
explaining BSs found in eccentric and/or wide binaries but 
is also required to increase the number produced. We can 
expect this theoretically through the hardening of primor- 
dial binaries, which increases the chance of mass transfer, as 
well as the possibility of collisions between MS stars in the 
cluster core. In addition, dynamical evolution together with 
the effects of a tidal field alters the cluster mass function 
as low-mass stars are stripped preferentially from the outer 
regions (Terlevich 1987). 

A problem with all of the population synthesis runs is 
that the number of RS CVn systems is comparable to the 
number of BSs, whereas the observations suggest that there 
should only be one for every four BSs. This is most likely an- 
other signature of the cluster dynamics. Giant stars present 
a much larger cross-section for collisions than MS stars and 
are therefore more likely to be involved in dynamical en- 
counters (even though their existence is shorter). Moreover, 
the hardening of close binaries and the disruption of wide 
binaries both act to reduce the number of RS CVn systems. 



ing initially, the number of BSs produced actually increases 
with time. This is because increasingly more stars are evolv- 
ing off the MS, growing in radius and interacting with their 
companion. Also, as time goes on, the masses of the pro- 
genitor stars that interact to produce BSs decreases so that 
the BS lifetime increases. After some time the number of 
BSs produced peaks. It then starts to fall, mainly as a re- 
sult of a decreasing number of progenitor systems. The time 
of this peak is largely dependent on the ratio of the mean 
binary separation to the size of a star at the MS turn-off. 
The latter decreases with time and if stars do not grow to 
such large radii then binary interaction is less likely. Addi- 
tionally, as this ratio increases, stars become more likely to 
fill their Roche-lobes on the GB, if at all, resulting in more 
cases of common-envelope evolution rather than steady mass 
transfer. We find the peak in TV B s for run PS6 occurs at 
about 5 Gyr quite unlike Pols & Marinus (1994) who found 
a peak at about 0.1 Gyr for a separation distribution flat 
in log a. Their mass function was biased towards stars with 
M > 2 Mq because they were interested in younger clusters. 
At 5 Gyr the turn-off mass drops below 1.25 Mq for the first 
time which means that the brighter MS stars begin to have 
radiative cores and therefore do not rejuvenate as much af- 
ter mass transfer. As a result, BS lifetimes decrease in re- 
lation to MS lifetimes, contributing to the decrease seen in 
TV BS /TV M s after 5 Gyr. Stars of lower metallicity have shorter 
MS lifetimes for M 9 Mq so that for a particular age the 
MS turn-off mass decreases with decreasing cluster metallic- 
ity. Run PS9 has a lower metallicity than PS6 and has a MS 
turn-off mass of 1.24 Mq at 4.2 Gyr. Consequently the peak 
in TV BS occurs at about this time which helps to explain why 
more BSs are produced at 4.2 Gyr than by PS6. Also, stars 
of the same mass grow to larger radii on the MS for lower 



4 THE TV-BODY CODE 

The study of star clusters is currently at a very exciting 
stage. Observationally the improved resolution of the Hub- 
ble Space Telescope (HST) is providing a wealth of high- 
quality information on clusters and their stellar populations 
(e.g. Guhathakurta et al. 1998; Piotto et al. 1999). This in- 
cludes the dynamically old globular clusters which populate 
our Galaxy as well as those of the Magellanic Clouds which 
exhibit a wide range of ages. Coupled to this are the recent 
advances in computer hardware which have brought the pos- 
sibility of direct globular cluster modelling within reach for 
the first time (Makino 1999). 

In practice the direct integration of an TV-body system 
presents many technical challenges and has only limited ap- 
plicability to real clusters because the required value of TV 
is too large for a simulation to be completed in a reasonable 
time. As a result, most TV-body simulations performed so 
far have involved a varying number of simplified and unre- 
alistic conditions, such as including only single stars, using 
only equal-mass stars, neglecting stellar evolution or assum- 
ing no external tidal field (e.g. McMillan, Hut & Makino 
1991; Heggie & Aarseth 1992). Despite the fact that, even 
with simplified conditions, direct TV-body simulations are 
limited to TV considerably less than is needed for globu- 
lar clusters, much can still be learnt by scaling the results 
with particle number, or by modelling small open clusters. 
The scaling of time depends essentially on the mechanism 
to be modelled (Meylan & Heggie 1997) but unfortunately 
the various timescales scale differently with TV so complica- 
tions arise when competing processes are involved. Aarseth 
& Heggie (1998) discuss these problems further and present 
a hybrid time-scaling method which can cope with the tran- 
sition from early evolution related to the crossing time to 
later relaxation-dominated evolution. 
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A complication in the use of the results of small- TV cal- 
culations to make inferences relating to larger clusters is that 
many of the structural properties are TV-dependent (Good- 
man 1987). Also energy considerations which may be dom- 
inant for small- TV, such as the binding energy of a single 
hard binary, may not be significant for a larger system where 
only the overall energetics is important. Another problem is 
that, as TV decreases, the results become increasingly noisy 
owing to statistical fluctuations. Casertano & Hut (1985) 
have studied how noise affects determination of the core pa- 
rameters while Giersz & Heggie (1994) have suggested that 
the situation may be alleviated by averaging the results of 
many simulations. The validity of the results of TV-body cal- 
culations has been challenged on a fundamental level by 
Miller (1964) who showed that two TV-body systems inte- 
grated from similar initial conditions diverge exponentially. 
Quinlan & Tremaine (1992) note that this instability occurs 
on a timescale comparable with the crossing time, mak- 
ing the results of TV-body integrations extremely sensitive 
to numerical errors and therefore unreliable over relaxation 
timescales. However, they also find that the numerical or- 
bits of TV-body models are shadowed by real systems even 
though their initial conditions differ. 

Ultimately the TV-body approach remains the method 
of choice for creating dynamical models of star clusters be- 
cause a minimum number of simplifying assumptions are re- 
quired and it is relatively easy to implement additional real- 
istic features. For a comprehensive description of the various 
methods for dynamical star cluster modelling see Meylan & 
Heggie (1997), or Hut et al. (1992). 

The TV-body code we use is NBDDY4, versions of which 
have been described in the past by Aarseth (1996, 1999a, 
1999b). This code has been adapted to run on the HARP-3 
special-purpose computer (Makino, Kokubo & Taiji 1993) 
and makes use of the many advances made in the field since 
publication of the original direct TV-body code (Aarseth 
1963). 

4.1 Integration 

In TV-body simulations it is usual to choose scaled units 
(Heggie & Mathieu 1986) such that G = 1, the mean mass 
is 1/TV and the virial radius is unity. This means that, in 
TV-body units, the initial energy of the system in virial equi- 
librium is —1/4 and the crossing time is 2\/2. The TV-body 
units can be scaled to physical units via the total stellar mass 
and an appropriately chosen length-scale factor. In this work 
we do not model the initial phase of cluster evolution, the 
formation of the cluster. Instead the initial conditions are 
chosen subject to suitable observational constraints and we 
then integrate a system in virial equilibrium forward in time. 
Basic integration of the equations of motion is performed by 
the Hermite scheme (Makino 1991) developed specifically 
for the HARP. This scheme employs a fourth-order force 
polynomial and exploits the fast evaluation of the force and 
its first time derivative by the HARP. It is more accurate 
than the traditional divided difference formulation (Ahmad 
& Cohen 1973; Aarseth 1985) of the same order and has 
advantages in simplicity and performance. 

To exploit the fact that stellar systems involve a wide 
range of particle densities, and thus that different parti- 
cles will have different timescales for significant changes to 



their orbital parameters, it is desirable to introduce indi- 
vidual time-steps. The individual time-step scheme requires 
the coordinate and velocity predictions to be performed for 
all other particles at each force evaluation. Fortunately the 
overheads that this introduces are reduced by utilizing the 
HARP hardware for fast predictions to first order. Signifi- 
cant gains in efficiency are also achieved by adopting quan- 
tized hierarchical time-steps. Although first developed by 
McMillan (1986) to optimize vectorization of TV-body codes, 
the method also aids efficient parallelization of the force cal- 
culations. The advantage of quantized time-steps is that a 
block of particles can be advanced at the same time so that 
only one prediction call to the HARP is required for each 
block. 



4.2 Close Encounters and Regularization 

When a binary system is within an environment such as a 
star cluster where gravitational encounters with other bod- 
ies are possible it is not sufficient to describe the orbit by 
averaged quantities, as used for isolated systems in binary 
population synthesis (Hurley, Pols & Tout 2000). Instead 
the orbit must be integrated directly, in a way that enables 
the positions of both stars to be known at any time. This 
is complicated by the fact that the orbital characteristics 
are affected by perturbations from the attractive forces of 
nearby stars. If the relative separation of the two binary 
stars is R then 

R=-^R + P (5) 

describes the motion where P represents the external tidal 
perturbations. As R — > this equation becomes strongly sin- 
gular, leading to increasing errors and small time-steps if in- 
tegrated by standard methods. Therefore alternative meth- 
ods which regularize the equations of motion for close en- 
counters, or make use of a hierarchical data structure, must 
be employed. This includes the case of two strongly inter- 
acting particles in a hyperbolic encounter. NB0DY4 makes use 
of KS regularization (Kustaanheimo & Stiefel 1965) which 
treats perturbed two-body motion in an accurate and effi- 
cient way. A recent development has seen the introduction 
of the Stumpff KS scheme (Mikkola & Aarseth 1998) which 
achieves a high accuracy without extra cost. This scheme, 
also incorporating the so-called slow-down principle based 
on adiabatic invariance, is now used in NB0DY4 and requires 
about 30 steps per orbit to obtain high accuracy when rel- 
atively weak perturbations are involved. 

Since all the KS integrations are carried out on the host 
machine rather than the HARP, and each orbit may require 
many KS steps, the treatment of many perturbed systems 
is often the most expensive part of the simulation. There- 
fore the extent of the primordial binary population, and its 
distribution of periods, is a major consideration when deter- 
mining the size of a simulation that can be completed in a 
reasonable time. For perturbed two-body motion, only con- 
tributions from relatively nearby particles need be consid- 
ered because the tidal force varies as 1/r 3 . If no perturbers 
are selected for a KS pair at apastron then unperturbed 
KS motion is assumed and the system can be advanced one 
or more periods without stepwise integration. Because hard 
(i.e. close) binaries are less likely to be perturbed, a signif- 
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icant number of such binaries among the primordial popu- 
lation can reduce the load on the host machine. Thus the 
choice of period distribution is once again a prime consider- 
ation. 

To cope with strong interactions between close bina- 
ries and single stars the KS treatment has been extended 
to three-body regularization (Aarseth & Zare 1974). This 
basically requires two KS regularizations coupled by suit- 
able coordinate transformations. A further generalization to 
include a fourth body led to the formulation of chain reg- 
ularization (Mikkola & Aarseth 1990, 1993). The essential 
feature of this challenging treatment is that dominant in- 
teractions along the chain of particles are modelled as per- 
turbed KS solutions and all other attractions are included 
as perturbations. This is ideal for treating compact subsys- 
tems and is currently used in NB0DY4 for configurations of 
three to six bodies. 

4.3 Stellar Evolution Treatment 

To provide realistic cluster models the simulation must be 
able to account for changes to the radii and masses of stars 
during the lifetime of a star cluster. The evolution of the 
single and binary stars must be performed in step with in- 
tegration of the dynamics so that interaction between the 
two processes is modelled consistently. The earlier version 
of NB0DY4 included stellar evolution in the form of the algo- 
rithms presented by Tout et al. (1997) but this treatment is 
only relevant to stars of Population I composition. It is also 
based on stellar models which have since been superseded by 
the models of Pols et al. (1998) incorporating, among oth- 
ers, improvements to the equation of state, updated opacity 
tables, and convective overshooting. We have now incorpo- 
rated the updated stellar evolution treatment described by 
Hurley, Pols & Tout (2000), valid for all metallicities in the 
range Z = 10~ 4 to 0.03, into NB0DY4 in its entirety. This 
means that the metallicity of our cluster models can be var- 
ied and that a more detailed and accurate treatment of all 
the single star evolution phases is used. Variations in com- 
position can affect the stellar evolution timescales as well as 
the appearance of the evolution in a CMD and the ultimate 
fate of a star. 

All stars carry with them a set of variables that describe 
their evolutionary state. These are the initial mass (Mo), 
current mass {Mt), stellar radius (R*) and stellar type (fc*). 
Other physical parameters, such as the luminosity and core 
mass, are not saved for each star because these are required 
less frequently and are calculated when needed. 

Because the stars evolve at different rates, depending 
on their evolutionary stage and mass, they require different 
frequencies for the updating of their variables. Each star 
has the associated variables T ev and T cv o which represent 
the next stellar evolution update time, and the time of the 
last update, respectively. Whenever T > T ov , where T is the 
simulation time, the star is updated by the stellar evolution 
algorithms. The value of T cv is then advanced by a suitable 
choice of time-step At which depends on the evolutionary 
stage of the star (see Hurley, Pols & Tout 2000) and ensures 
that the physical parameters change sufficiently smoothly. 
After At is chosen we check whether the stellar radius will 
change by more than 10% over the interval and reduce At if 
this is the case. In practice the stellar evolution algorithms 



are not called continuously because successive calls to the 
main routine which controls the treatment are limited by 
a minimum time interval, which is fairly small. This means 
that one call to the routine can involve the updating of many 
stars and that an individual star may be advanced several 
times within one call if it is in a particularly rapid stage of 
evolution. 

In the simulation all stars have evolved for the same 
amount of time, i.e. the age of the cluster, but they may 
have different relative ages, for example if a star has been 
rejuvenated by mass transfer. Also, when a remnant stage 
is begun, such as a helium star, white dwarf (WD), neutron 
star (NS) or black hole (BH), the age of the star is reset 
to zero for the new type. Therefore the additional variable 
EPOCH is introduced for each star so that its age at any time 
is given by the difference between the current cluster age and 
its EPOCH. 

We include mass loss by stellar winds according to the 
prescription given in Hurley, Pols & Tout (2000). The time- 
step is limited so that a maximum of 2% of the stellar mass 
can be lost. Any mass lost which is not accreted by a close 
companion is assumed to leave the cluster instantaneously, 
with the appropriate corrections made to account for the 
change in potential energy, maintaining energy conservation 
for the system. The force on nearby stars must also be modi- 
fied as a result of the mass change or, if a significant amount 
of mass is lost, a complete re-initialization of the force poly- 
nomials on the HARP may be required. Also, if a NS or 
BH is formed, a velocity kick taken from a Maxwellian dis- 
tribution with dispersion a = 190 km s -1 (see Hurley, Tout 
& Pols 2000) is given to the supernova remnant which is 
usually enough to eject it from the cluster. 

4.4 Binary Evolution Treatment 

For the treatment of evolution within a binary system we in- 
clude the features described in Hurley, Tout & Pols (2000) 
over and above those of the old algorithms (Tout et al. 1997). 
However, we use the treatment of tidal circularization devel- 
oped by Mardling & Aarseth (2000, hereinafter MA) which 
includes features that cope with the added complication of 
external perturbations to the orbital parameters, and has 
been implemented to work closely with the KS scheme. 

Each binary has an associated composite particle, the 
centre-of-mass (CM) particle, which has its own set of vari- 
ables, Mt, k», T cv etc., as well as additional variables such 
as the binding energy per unit mass of the system. A typical 
binary will make a journey through a series of CM evolution 
stages, beginning with fc»,cM = when it is first created: 
either primordially or during the evolution. A close binary 
orbit circularizes tidally as energy is dissipated but angular 
momentum is conserved. Tidal circularization (fc*,cM = —2) 
is activated when the circularization timescale of the binary 
is less than 2 x 10 9 yr. Note that the intrinsic stellar spin is 
not incorporated in the circularization treatment so that the 
angular momentum of an unperturbed binary orbit remains 
constant during the circularization process. This could be 
rectified in line with the tidal model of Hurley, Tout & 
Pols (2000). 

For a binary of high eccentricity it is possible for the en- 
ergy exchange between the orbit and the tides at periastron 
to become chaotic (k*,cM = — !)• This means that oscilla- 
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tions occur which gradually damp the system until sufficient 
energy is dissipated for the chaos boundary to be crossed, 
when the orbit settles on some eccentricity whence angular 
momentum is conserved, and the orbit begins to circularize 
(MA). The implementation is based on a non-linear dissi- 
pation timescale and assumes conservation of total angular 
momentum for the system (i.e. the stars spin up). Chaotic 
behaviour is most likely to occur in tidal-capture binaries 
and MA note that it is still not clear how stars in chaotic 
orbits respond to the huge tidal energies involved. Tout & 
Kembhavi (1993) and Podsiadlowski (1996) have shown that 
the response of the tidally heated primary depends on the 
region of the star in which the energy is deposited and also 
the timescale on which this energy is thermalized within the 
star. 

Once the circularization is complete (fc*,cM = 10) the 
binary is tested regularly for Roche-lobe overflow (RLOF). 
Prior to circularization the unlikely possibility of RLOF is 
ignored. In most cases this technicality does not lead to 
any physical irregularities because tidal interaction gener- 
ally acts to remove any eccentricity on a timescale shorter 
than the evolution timescale of the binary. However prob- 
lems may arise if a system forms in a close eccentric orbit. 
Even then there is evidence that some circularization occurs 
during the formation process (Mathieu 1994). While the or- 
bit is circularizing the stars are also evolving and growing 
in radius, so contact is expected in some systems as the or- 
bit approaches zero eccentricity. However the timescale for 
a collision is actually prolonged because the periastron dis- 
tance grows as e shrinks, owing to conservation of angular 
momentum. Thus, as the stars are expanding so is the min- 
imum separation between them. Also most binaries do not 
begin to circularize, or to come close to a RLOF state, until 
one of the stars is on the GB or AGB, at which point it is 
losing mass in a stellar wind that causes an additional in- 
crease in the separation. Nevertheless, stages of rapid radius 
growth may result in premature coalescence. 

It is also possible for the eccentricity of a binary to in- 
crease as a result of external perturbations acting on the 
orbit. If this induces an eccentricity in an orbit which was 
previously circular the binary is reset to standard type 
(&*,cm = 0). This is also done if an eccentric binary sur- 
vives one of the component stars exploding as a supernova 
to leave a NS or BH remnant. 

While a binary is in a detached state the stars are up- 
dated in the same way as single stars. The only difference 
is that both stars must be treated at the same time so that 
the possibility of wind accretion can be modelled. If the CM 
type is fc*, C M = 10 or greater, each time the binary stars 
are updated the time until the primary star fills its Roche- 
lobe is estimated. This is used to set T eViC M- If at any stage 
T > T ev ,cM then RLOF could have begun and the binary 
is subject to the roche procedure (fc*,cM = 11) where it 
is evolved forward according to the treatment described in 
Hurley, Tout & Pols (2000). This means that the physical 
time of the binary evolution can move slightly ahead of the 
cluster integration time, introducing the need for coasting 
periods in which the binary is put on hold until the dy- 
namical integration time catches up. The amount of time 
that the binary moves ahead during ROCHE is determined 
by physical conditions of the components. Any stellar wind 
mass loss during RLOF is also dealt with inside ROCHE so 



it is important that T cv for each of the component stars is 
greater than the time reached by the binary when it exits 
the active RLOF stage, i.e. T C v,cm < Tev,i,T eVi 2- In this 
way the component stars have their evolution treated as 
part of the RLOF process and not as individuals. A sus- 
tained phase of RLOF may involve a series of active ROCHE 
calls each followed by a coasting period. It is also possi- 
ble that a binary may experience more than one phase of 
RLOF during its evolution. This is easily accommodated by 
the treatment. During mass transfer the separation of the 
binary stars changes and the KS variables are updated fre- 
quently. As a fundamental variable for two-body regulariza- 
tion the binding energy per unit mass must be determined 
accurately at all times. In fact, all the KS variables must be 
corrected for mass loss. A mass-loss correction is also made 
to the total energy of the cluster and the force polynomials 
of nearby stars are re-initialized. 

We treat common-envelope evolution and collisions that 
arise as a result of contact systems as part of the ROCHE pro- 
cess. In general, direct stellar collisions during the cluster 
evolution are very rare because the interacting stars most 
likely form a tidal capture binary, which may be quickly 
followed by coalescence via a common-envelope or contact 
phase. This is due to the relatively low velocity dispersion, 
a ~ lOkms - , of cluster stars. In galactic nuclei where 
typically a ~ 100 kms -1 direct collisions are expected. We 
use a periastron criterion, derived for main-sequence stars 
(Kochanek 1992), 

where Ri is the primary radius, to determine direct colli- 
sions. 

Because common-envelope events and collisions fre- 
quently lead to coalescence, either one or two particles can 
become redundant. Similarly some single star supernovae 
may not leave a remnant. In such cases it is simplest to 
create a massless component and place it well outside the 
cluster so that it escapes soon. If the component of a binary 
is removed then the coalescence product is given the CM 
coordinates and velocity of the original binary. Energy cor- 
rections associated with mass loss are performed and a new 
force polynomial initialized. 

As already mentioned, the fraction of the stars in pri- 
mordial binaries is of great importance. It is generally ac- 
cepted that some degree of primordial population is present 
because the rate of binary formation during the evolution 
would not be enough to halt core-collapse appreciably (Hut 
et al. 1992). Binary stars in clusters can also be identified 
through their position in the CMD of the cluster as the 
combination of light from the two stars displaces the binary 
from the position of a single star having the same mass as 
either of the components (Hurley & Tout 1998). This is most 
noticeable on the MS where the existence of a distinct bi- 
nary sequence has been exploited by the resolution of the 
Hubble Space Telescope (HST) to reveal globular cluster bi- 
nary fractions of 10 to 30% (e.g. Richer et al. 1997; Elson 
et al. 1998). For open clusters results are often uncertain 
because the number of stars is smaller and membership can 
be difficult to determine. Even so, a significant number of 
pre-main-sequence binaries and multiple systems have been 
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found in young open clusters and associations (e.g. Simon 
et al. 1995; Brandner et al. 1996). 



4.5 Hierarchical Systems 

During the evolution of a star cluster stable multiple sys- 
tems can form as a result of dynamical interactions. The 
formation of hierarchical triples, in which one component of 
a binary is itself a binary, has been shown to occur in var- 
ious scattering experiments (e.g. Mikkola 1983; Bacon, Sig- 
urdsson & Davies 1996) mainly as a result of strong binary- 
binary encounters. The inner binary of such a system may be 
relatively hard so that integration of the strongly perturbed 
KS solution can prove to be extremely time-consuming. This 
is especially true if the hierarchy is stable for a long pe- 
riod of time. In that case the inner binary only experiences 
short-term fluctuations in its orbital parameters, and so it 
is acceptable to perform direct integration only of the outer 
orbit while the system remains stable. 

A semi-analytical stability criterion based on the anal- 
ogy with the chaos boundary in tidal evolution has been 
presented by MA and implemented in NB0DY4. They employ 
a critical periastron distance for the outer orbit in terms of 
the inner semi-major axis, ai n , given by 
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where e out is the eccentricity of the outer orbit and C ~ 2.8 
is determined empirically. The mass-ratio of the outer or- 
bit is g out = Ma/ {Mi +M 2 ), where Mi and M 2 are the 
component masses of the inner orbit. If the outer periastron 
separation is greater than -Rp"*rit we consider the hierarchy 
to be stable and temporarily merge the stars into one KS 
system consisting of the CM particle of the inner binary and 
the third body, M3. Because the period of the outer orbit is 
considerably longer than that of the inner orbit this proce- 
dure greatly reduces the computational cost. Quadruple and 
higher-order systems are similarly dealt with. Procedures to 
model the cyclic oscillations of the inner eccentricity, the 
Kozai effect (Kozai 1962), and any tidal circularization that 
is induced, are implemented according to MA. 

The possibility of an exchange interaction, in which one 
of the inner binary components is displaced by an incoming 
third star, is checked by the criterion of Zare (1977)0. As 
noted by Aarseth (1999a), the stability boundary lies above 
the exchange boundary when all the masses involved are 
comparable and the two only begin to overlap when g out — 5. 
If an exchange does occur then the expelled star invariably 
leaves the three-body system altogether. 

To ease the book-keeping required when particles are 
regularized or hierarchies are formed, the simulation parti- 
cles are kept in an ordered data list throughout the evolu- 
tion. Consider a model composed of iV s single stars and iVb 
regularized binaries, i.e. N — N s + 2Nb stars in total. The 
first 2iVb entries in the data list are the binary stars, with 
each binary pair grouped together, followed by the N s single 
star entries, and finally the Nb CM particles. Take the case 



t Because of the degeneracy of angular momentum, this criterion 
is only used for small inclinations 



of a basic KS binary which is combined with a single star 
to form a stable hierarchy. The single star is moved to the 
position formerly occupied by the second component of the 
binary and the CM particle of the old binary is moved to 
the first component position. The values, such as binding 
energy, of the old binary CM are moved to the position for- 
merly occupied by the single star which now has zero mass 
and is a ghost particle, i.e. it does not contribute to force cal- 
culations. Component masses of the inner binary are saved 
in a merger table and the CM position for the binary now 
holds variables relevant to the outer orbit. If the hierarchy 
is broken up then all the variables are easily re-assigned and 
the inner binary restored. Possible reasons for termination 
of a hierarchy include violation of the stability criterion or 
the onset of RLOF in the inner binary. If the hierarchy in- 
volves the merger of two binaries then both CM particles 
would be the components of the new binary. 



5 ESCAPE FROM A TIDALLY-LIMITED 
CLUSTER 

Fan et al. (1996) present CCD spectrophotometry of 6 558 
stars in the field of M67. They estimate that the MS is com- 
plete down to 0.5Mq and that the cluster has a total ob- 
served mass of approximately 1 OOOM0 in stars with masses 
greater than 0.5Mq. The observations reveal a mean (pro- 
jected) half-mass radius of 2.5 pc for MS stars compared 
with 1.6 pc for the (more centrally condensed) BS popula- 
tion. The tidal radius is at 10 pc and their data is consistent 
with a 50% binary fraction. 

With an initial population comprising single star masses 
chosen from the KTG3 IMF, binary masses from the KTG1 
IMF and a 50% binary fraction, about 7000 stars (3 500M©) 
are required to give a current cluster mass of 2 5OOM0 (corre- 
sponding to 1 000M Q observed above 0.5M Q ). The mass loss 
is due solely to stellar evolution and assumes no dynamical 
effects on the population. However, in a cluster environment 
stars undergo gravitational interactions so that from time to 
time an encounter gives enough energy to a star that it can 
escape from the system. The timescale over which the stars 
evaporate in this way is related to the relaxation timescale, 
t t , of the cluster and is shortened by the presence of an ex- 
ternal tidal field as well as the inclusion of a full spectrum 
of stellar masses. 

An isolated spherical cluster with potential $ has an 
escape speed v e at radius r given by 



-2$ (r) 



(8) 



The mean-square escape speed in a system with uniform 
density is 



<^ 2 > = 4(^ 2 ) , 



(9) 



so the RMS escape speed is twice the particle RMS speed. 
For a Maxwellian velocity distribution the fraction of par- 
ticles that have speeds exceeding twice the RMS speed is 
7 = 7.4 x 10" 3 (see Binney & Tremaine 1987, p. 490) and 
the evaporation process can be approximated by removing 
a fraction 7 of stars each relaxation time, i.e. 



dN 
dt 



N 



(10) 
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Evaporation then sets an upper limit to the lifetime of the 
bound stellar system of about 10 2 i r . Note that numerical 
solution of the Fokker-Planck equation gives 7 = 8.5 x 1CF 3 
(Spitzer 1987, p. 54). The total energy of a cluster with total 
mass M and radius R is, according to the virial theorem, 



E=-k 



GM 
R 



(11) 



where k is a dimensionless constant of order unity. If it is 
assumed that the cluster evolution is self-similar, so that its 
shape remains fixed, k is independent of time. Since evapo- 
ration is mostly driven by weak encounters the stars escape 
with very little energy so that E remains essentially fixed. 
Therefore 

r Q \ M 1 

and the cluster contracts as it loses mass. 

The evolution of a real star cluster is somewhat dif- 
ferent to that of an isolated uniform system because the 
cluster is subject to the tidal force of the galaxy in which 
it resides. Consider a cluster in a circular orbit around a 
galaxy at a distance R G from the galactic centre. To esti- 
mate the strength of the tidal field exerted on the cluster 
we assume that M G , the galactic mass enclosed by the or- 
bit, is distributed throughout the inner spherical volume. 
By choosing the origin of a rotating reference frame to be 
the cluster centre-of-mass, the a;-axis directed away from the 
galactic centre, the j/-axis in the direction of rotation, and 
linearizing the tidal field, the equations of motion are 



2uj G y - 

2uj g x 
p 2 



o 2 
6UJ r X 



(Giersz & Heggie 1997), where 



Ug = 



GM G 
R G 



(13) 
(14) 



(15) 



is the angular velocity. The tidal radius r t of the cluster is 
defined by the saddle point on the z-axis of the effective 
cluster potential, analogous to the definition of the Roche- 
lobe radius in a binary system (except that the cluster is 
not co-rotating with its orbit). 

The Galactic tidal field can conveniently be described 
in terms of Oort's constants, 

A = 14.5 ± 1.5 kms^kpc" 1 
B = -12±3kms _1 kpc _1 

(Binney & Tremaine 1987, p. 14), so that 

A-B = (^-) = 26.5±4kms~ 1 kpc~ 1 . (16) 

This is consistent with an orbital velocity of v G = 220km s _1 
(Chernoff & Weinberg 1990) at R = 8.5 kpc. In this formu- 
lation the tidal radius is 



n = 



GM 



1/3 



dM 
~dt 



4A (A - B) 

The cluster escape rate can be expressed by 
k c M 



log in A t rh 



(17) 



(18) 



. _ M 4925 M 
M Q 4970 M 
M n 476 M, 



no binaries r h Q 2.34 pc 

1000 binaries r h Q 2.32 pc _ 

1000 binaries r h Q 1.39 pc 

330 binaries r h Q 0.22 pc_ 



tArh 

Figure 5. Evolution of the total mass of bound stars for TV-body 
models started with 10 000 stars and a standard tidal field with 
the time scaled by the current half-mass relaxation time of the 
cluster. Also shown is the evolution for a model started with only 
1 000 stars but with a 50% binary fraction. 



where M is the cluster mass at time t for N stars, A = QAN 
is used in the Coulomb logarithm and 



t A = 0.894 



N 



3/2 



log 10 A MV2 



Myr 



(19) 



is the half-mass relaxation time if rh is in pc and M in Mq. 
The constant k e quantifies the rate of mass lost in stars 
stripped from the cluster, where 7 = fc c /log 10 A because 
M oc N. Numerical solution of the Fokker-Planck equation 
for a tidally truncated cluster gives 7 = 4.5 x 10 -2 (Spitzer 
1987, p. 59). 

We have investigated cluster escape rates for simula- 
tions with a full mass spectrum, mass loss from stellar evo- 
lution and a standard Galactic tidal field (v G — 220kms~ 
and R G = 8.5 kpc), in a series of A^-body models evolved 
with NB0DY4. In these simulations we assume the stars are 
stripped from the cluster when r > 2r t . This underesti- 
mates slightly the escape rate compared to letting them 
escape when they are outside the tidal radius (Giersz & Heg- 
gie 1997). We performed a total of six simulations starting 
with N = 10 000. In each model the initial positions and 
velocities of the stars are assigned according to a Plummer 
model (Aarseth, Henon & Wielen 1974) in virial equilib- 
rium with the tidal radius defined by eq. (^)- The single 
star masses are chosen from the KTG3 IMF and the binary 
masses from the KTG1 IMF, with minimum and maximum 
single star limits of 0.1 and 50 Mq . We set the distribution of 
mass-ratios for the binaries to be uniform between 0.1 and 
1 and the metallicity at Z — 0.02. We choose the binary 
separations from the same distribution as population syn- 
thesis run PS1. Three of the simulations had no primordial 
binaries and a length scale chosen so that rh,o — 0.1rt,o- The 
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M n 4800 



no binaries r, Q 2.34 pc 



M 4925 M Q 1000 binaries r h 2.32 pc 
M Q 4970 M Q 1000 binaries r h 1.39 pc 



0.4 



0.2 



M/M„ 



Figure 6. Evolution of the maximum and half-mass radii for the 
same N = 10 000 models as in Figure |^. Note that the maximum 
radius is given by the radius of the star closest to the tidal radius. 
It quickly approaches the tidal radius following the initial expan- 
sion of the cluster. Since stars are not removed from the cluster 
until r > 2rt there is a small population of cluster members with 
7"max < r < 2rv On average the mass in cluster stars outside the 
tidal radius is less than 2% of the total cluster mass. 

other three had a primordial binary fraction of /b = 0.1 
two of these with the same initial half-mass radius as the 
single star models and one, started with rh,o — 0.06r t , more 
condensed. In each simulation the initial cluster mass, Mo, 
is about 5OOOM0 which gives a tidal radius of about 22 pc. 

The evolution of the cluster mass as a function of the 
number of half-mass relaxation times is shown in Figure || 
The inclusion of binaries has little effect, nor does varying 
the initial concentration of the cluster, as long as the time 
is scaled by the current half-mass relaxation time. In fact 
the initial size has little effect because, owing to mass loss 
from massive stars and core contraction, the cluster quickly 
expands to fill its tidal radius which depends only on mass. 
This is demonstrated in Figure ^ which shows the varia- 
tion of the half-mass radius with cluster mass as well as the 
position of the star closest to the tidal radius at the time. 
As each cluster evolves, the half-mass radius, which initially 
expands freely, starts to feel the effect of the contracting 
tidal radius so that the rate of expansion slows. Eventually 
the half-mass radius begins to contract and the cluster evo- 
lution becomes remarkably self-similar. For each simulation 
the turn-over occurs at M/Mo ~ 0.5 when t/i r h — 6. During 
the subsequent self-similar evolution phase the ratio of the 
tidal and half- mass radii is ft /rh — 4. 

The mean core-collapse (CC) time for all the models 

t We define the binary fraction as /b = iVb/(iVb + N s ) for JVb 
binaries and N B single stars. During this work we may also refer 
to a binary fraction as a percentage, i.e. /b = 50%, and we note 
that this should more correctly be called a binary frequency. 



is t C c — 4t r h. In terms of the initial half- mass relaxation 
time this is tec — 7i r h,o because t r h,cc — 2.2t r h,o owing to 
the initial expansion of r^. At the time of core-collapse an 
average of five hard binaries have formed in the core of the 
single star models. The model with the smaller rh,o reached 
the rh turn-over point in about 85% of the time it took the 
lower density models but this single result is not significant. 

From the N = 10 000 models we find that k c ~ 0.3 fits 
the data adequately for almost the entire cluster evolution. 
This result is consistent with a series of six N = 1 000 sim- 
ulations that we began with /b = 0.5 and r^o — 0.03r t . For 
these smaller models the half-mass radius starts to contract 
at M/Mo — 0.5. The mass-loss rate is also compatible with 
the results of Giersz & Heggie (1997) from a series of iV-body 
models with N = 500, when we take into account the fact 
that Plummer models have a lower escape rate than King 
models. This is due to Plummer models having a slightly 
weaker tidal field, and therefore evolving more slowly in the 
post-core-collapse phase of evolution, than their King model 
counterparts. Giersz & Heggie (1997) use a minimum mass 
of 0.4 Mq, compared to the value of 0.1 Mq in the models 
presented here, so they have a higher proportion of massive 
stars. Their shorter evolution timescales cause the simula- 
tions to evolve faster which increases the escape rate and 
partly explains why their value of k c ~ 0.6 (converted to the 
units used here assuming A = 0.4A for both sets of models) 
is larger. Giersz & Heggie (1997) found that A = 0.015JV 
is required to make the results of their models agree with 
the Fokker-Planck models of Chernoff & Weinberg (1990), 
in which case the conversion of their escape rate gives an 
even larger value of k c ~ 1.5. 

So how robust is our k c against changes to the initial 
conditions assumed for the Af-body models? We have al- 
ready discussed that it is sensitive to the model from which 
the density and velocity profiles are taken: Plummer, King 
or otherwise. The mass function used also has an effect even 
though a change in IMF slope will be offset to some de- 
gree by the inverse dependence of the relaxation timescale 
on average stellar mass. Our escape rate has been derived 
by considering a number of model sizes and binary fractions 
but must be tested over a greater range of parameter space 
before it can be universally accepted. Factors such as the 
proportion of hard binaries must also have an, as yet unde- 
termined, effect on the result. Caution should be exercised 
when applying k e to situations where the initial conditions 
differ substantially with those from which it was derived. 

If the evolution of the cluster is self-similar then 



rh,o 



\MoJ 



2-C 



(20) 



which stems from writing the rate of change of the total 
cluster energy as 



dE 
dt 



Cg dM 

U~dT 



(21) 



If £ > 2, rh increases as M decreases. This can occur during 
the initial violent relaxation phase when stellar wind mass- 
loss causes an overall expansion, or during the post-core- 
collapse expansion of the inner regions. If f = 5/3 then 
the evolution of the half-mass radius follows the decrease in 
tidal radius as the cluster evolves. Using eq. ( f?(i| ) to integrate 
cq- © gives 
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M (t) = Mo 



l - 



7 - 3C fee 



lo gio A o irh.O 



(22) 



We can use this with k e — 0.3 to estimate the initial mass 
and half-mass radius required to give the current values of 
M ~ 2 5OOM and r h ~ 2.5 pc at t ~ 4 200Myr for M67. 
However, as Figure ^ shows, a complication arises because 
£ as it is defined in eq. ( po) ) is not constant throughout 
the cluster lifetime. Furthermore, the evolution of a tidally- 
limited cluster is not self-similar initially. 

The initial cluster values can be roughly constrained 
by the following method. For M/M < 0.5 the TV = 10 000 
models show that f = 5/3 and r t = 4rv So we choose a 
value for Mo and use eq. ( |r?| ) to find r t i/ 2 , i.e. tidal radius 
when M = 0.5Mo. This defines rj, 1/2- Next we choose a 
value for rh,o which can be used with r^i/a in eq. (|o|) to 
calculate an approximate £ for the 1 > M/Mq > 0.5 phase of 
evolution. Eq. (|ltj) with TVo — 2Mo / Mq (assuming the av- 
erage stellar mass at t = 0.0 is roughly 0.5 Mq , which is true 
for the KTG3 IMF) gives t r h,o which we use in eq. (^2|) to 
calculate <a/2> the time taken for the cluster mass to reduce 
to half its initial value. From then on the initial parameters 
in eqs. (pc| ) and (53) can be replaced by the corresponding 
values at t]/ 2 to find M and as a function of time with 
£ = 5/3. By iterating on this method we can find the initial 
values corresponding to current observed cluster properties. 
It should be noted that the solution is not single- valued: 
increasing Mo and decreasing rh,o gives similar values. 

There is a problem with this method for M67, demon- 
strated if we put M = 2 500 M in eq. (0) for the stan- 
dard Galactic tidal field. This gives a tidal radius of 17.5 pc, 
corresponding to rh = 4.4 pc which is almost twice what 
is observed. Francic (1989) estimates a lower limit of 9pc 
for the tidal radius of M67 and finds no stars with high 
membership probabilities outside this range. Additionally 
the data of Fan et al. (1996), from which the current mass 
is estimated, extend no further than 10 pc from the cluster 
centre. This agrees well with the observed half-mass radius 
of 2.5 pc and rt = 4rh from the model data. So does M67 
contain less mass than is observed or does its tidal radius 
not correspond to the standard Galactic tidal field? Cher- 
noff & Weinberg (1990) find that taking v G = 220 km s^ 1 
for 3 < Rg < 20kpc is consistent with current theoretical 
models of the mass distribution of the Galaxy. The posi- 
tion of M67 relative to the Sun gives R G — 9kpc so the 
local tidal field should apply. However, the Galactic orbit of 
M67 is slightly eccentric (Carraro & Chiosi 1994), with an 
apogalacticon of 9.09 kpc and a perigalacticon of 6.83 kpc, so 
it has been subject to a time varying tidal field. Possibly the 
structure of M67 was altered by an event in its past, such as 
an interaction with another cluster or an interstellar cloud 
(Terlevich 1987). We do not dwell on this non-standard tidal 
radius, preferring to discuss it further in Section What is 
important is that the mass assumed for M67 corresponds 
to the observations from which it is derived, i.e. 2 500 Mq 
within 10 pc. 

Using this method we estimate that M67 had TVo — 
40 000 and rh,o — 1 pc to evolve to its current observed 
parameters. This is unfortunate because a simulation with 
No > 20 000 takes a prohibitively long time with currently 
available equipment, especially when using a 50% binary 
fraction. On the other hand the main interest of this work 
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a 0.06 in one case 

Table 3. Model parameters for all simulations described in this 
work. The final row lists the number of times each simulation was 
performed. All simulations share the following general properties 
(see text for details): (i) single star masses taken from KTG3 IMF 
with minimum and maximum mass limits of 0.1 and 50 Mq; (ii) 
binary masses taken from KTG1 IMF with a uniform distribution 
in q; (iii) thermal eccentricity distribution; (iv) Z = 0.02; and (v) 
escape radius of 2rv 



is TV B s near 4 200 Myr so it should be possible to extract 
meaningful results using a semi-direct method. 



6 M67 TV-BODY MODEL 

We start the simulation after the self-similar phase has be- 
gun, using an initial model that reflects the preceding stel- 
lar, binary and cluster evolution. The model is substantially 
smaller in population number than the TV ~ 40 000 required 
to start from the birth of the cluster, so it can be evolved 
directly to the age of M67 in a reasonable time. In view of 
the discussion of Section ^| and the small number statistics 
of the BS population, this is preferable to scaling up the re- 
sults of a smaller simulation starting from initial conditions. 
For a first attempt using this method we begin a simulation 
with TV = 15 000 stars at 2 500 Myr. The initial parameters 
for this semi-direct M67 model and for the simulations from 
which the escape rates were derived (see previous Section) 
are summarized in Table ^. 

To estimate the binary fraction and the maximum or- 
bital separation for the binaries in this starting model we 
consider the hard binary limit for the cluster. If indeed 
TVo = 40 000 and r h ,o = lpc then a hard ~ 10 AU for M67 
initially, whereas it is now more like 140 AU. Only binaries 
which were initially hard could survive the early phases of 
cluster evolution so a conservative estimate of 50 AU is used 
for the upper limit of the separation distribution in the start- 
ing model. Lowering this limit would increase the incidence 
of interaction within binary systems and presumably lead 
to a larger number of BSs produced by the simulation. This 
would not be accompanied by an increase in the number of 
BSs found in wide eccentric orbits as the population of wide 
binaries available for exchange interactions is decreased. As 
already mentioned, observations of M67 indicate that the 
current binary fraction could be as high as 50%. During the 
cluster evolution soft binaries are broken-up in dynamical 
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mass in 


mass in 


mass 
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fraction 


single stars 


binaries 


0.100 


0.137 


0.081 
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0.0 
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0.263 
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0.363 
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236.9 


0.503 


0.694 


0.193 


362.3 


785.6 


0.695 


0.960 


0.148 


378.3 


827.6 


0.961 


1.326 


0.116 


272.1 


1038.1 


1.327 


1.833 


0.055 


118.6 


732.5 


1.834 


2.532 


0.005 


16.9 


83.6 


2.533 


3.500 


0.004 


0.0 


98.6 



MF 

2.5 Gyr evolved MF 
2.5 Gyr dynamical MF 



Table 4. Mass groups used to populate the M67 starting model. 
All masses are in units of Mq. The lower and upper mass limits 
of each group are given in the first two columns, followed by the 
number fraction of stars in that group for the combined N = 
10 000 with / b = 0.11 models at 2 500Myr. 5 000 single stars 
and 5 000 binaries, all evolved to an age of 2 500Myr (see text 
for details), are chosen for the M67 starting model according to 
these number fractions. The resulting masses in single stars and 
in binaries for each mass group in the M67 starting model are 
given in the final two columns. 



encounters so this is a lower limit on the primordial frac- 
tion. However, hard binaries are retained preferentially by 
the cluster because they have a higher average mass than 
single stars. We therefore assume that the binary fraction 
has remained roughly constant during the cluster lifetime, 
in good agreement with the N = 10 000 data, so we use 
/b = 0.5 for the starting model. 

We use results from the previous simulations to estimate 
what the mass function (MF) will look like at 2 500 Myr. 
Figure @ shows the IMF for the N = 10 000 models that 
included binaries and the corresponding MF for the popula- 
tion at 2 500 Myr according to the population synthesis code 
and from the TV-body simulations at the same age. It is evi- 
dent that dynamics and the tidal field alter the MF, lowering 
it at the low-mass end and increasing the relative number 
of more massive stars. The number of BSs produced will be 
sensitive to the shape of the MF assumed for the starting 
model as a larger proportion of systems with mass compara- 
ble to that of the MS turn-off at a particular time will lead 
to more stragglers at that time. To generate the starting 
population for the semi-direct M67 simulation, we evolve a 
large population with a 50% binary fraction and Z = 0.02 
to an age of 2 500 Myr using the population synthesis code 
alone. The single star masses are chosen from the KTG3 
IMF and the binary masses and parameters are chosen in 
the same way as for the population synthesis run PS6, but 
with an upper limit of 50 AU in the separation distribution. 
Then using the information gained from the dynamically al- 
tered MF in Figure [j], i.e. the ratio of stars in each mass bin 
between the dynamical and non-dynamical MFs, we take 
5 000 single stars and 5 000 binaries from the large popula- 
tion to populate the starting model (see Table ^). This gives 
a starting mass of 6OOOM0 for the cluster at 2 500 Myr. 

We use a tidal field with v G = 350kms _1 at R G — 
8.5 kpc which fixes r t = 17 pc. This tidal radius determines 
the length scale used in the A-body simulation. The stars 




Figure 7. The full line shows the normalized IMF for a popu- 
lation with 11% binaries where the single star masses are chosen 
from the KTG3 IMF and the binary masses from the KTG1 IMF. 
The dashed line shows the mass function (MF) of the same popu- 
lation evolved to 2 500 Myr (note this is hidden by the full line for 
logm/Mo < —0.35) and the dash-dot line shows the MF of the 
same population evolved to 2 500 Myr in the iV-body code with a 
standard tidal field. 



M67 t = 2500.0 Myr 

N = 10000 t/t rh = 0.0 

N = 10000 t/t rh = 3.0 

N = 10000 t/t rh = 6.0 

N = 10000 t/t rh = 20.0 




Figure 8. Profile of the average stellar mass in successive La- 
grangian shells. The radius is scaled by the cluster half-mass ra- 
dius. Each Lagrangian shell contains 10% of the cluster mass. The 
profiles plotted are for the combined N = 10 000 models data at 
t/tih = 0.0, 3.0 (near core-collapse), 6.0 (when the model has lost 
half its mass) and 20.0, as well as for the initial M67 model at 
2 500 Myr. 
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are distributed according to a multi-mass King model (King 
1966; Chernoff & Weinberg 1990) with W = 7 and a central 
number density no = 2.28 x 10 3 pc -3 . The concentration of 
the model is determined by the dimensionless parameter 



Wo 



rr2 



(23) 



where is the central potential and we use a central ve- 
locity dispersion a 2 = 3 x 10 10 cm 2 s -2 . These positions are 
scaled so that the cluster just fills the tidal radius, which 
results in a half-mass radius of 3.4 pc for the starting model. 
Figure [] shows the evolution of the average stellar mass pro- 
file for the combined N — 10 000 model data. It can be seen 
that as the models evolve, two-body effects cause the heav- 
ier stars to segregate towards the inner regions. Using the 
King model to determine the initial spatial distribution of 
the stars for the M67 simulation builds in a degree of mass- 
segregation so that this energy equipartition is taken into 
account. 

We can estimate the relation between the mass of our 
starting model and its mass at 4 200 Myr using the escape 
rates discussed in Section ^, but this depends on how bina- 
ries contribute to the value of TV used in the calculation. If 
we assume that relatively hard binaries act as single stars 
when modelling relaxation effects then N would be 10 000 
and the mass at 4 200 Myr should be about 1 500 M . On 
the other hand, if the binary components behave as single 
stars then N = 15 000 and the mass left would be about 
3 000 Mo. Either way this is close enough to the mass de- 
rived from observations to proceed with the simulation. 

We evolve the model to T = 4 310 Myr using NB0DY4. 
At this time the cluster mass is 1 140 M©. It consists of 560 
single stars and 740 binaries, the tidal radius is 10 pc and 
the half-mass radius is 2.5 pc. The mass in single stars and 
binaries with masses greater than 0.5 M© is 1 050 Mq, more 
than 90% of the total mass. The simulation lasted for 1 260 
iV-body time units and took one month dedicated use of the 
HARP-3. 

Figure ^| shows the evolution of the number of BSs 
present in the simulation. Also shown are the numbers of 
triple systems, RS CVn systems and binaries containing a 
BS. At T = 4 200 Myr there are 22 BSs in the model but 
only one of these is in a binary. There is only one RS CVn 
system in the cluster at this time. The highest number of BSs 
present at any one time is 29 at T = 3 653 Myr with seven 
of these in binaries. At this time there are three RS CVn 
systems. 

Shown in Figure [l(] is the period-eccentricity distribu- 
tion of all BS binaries formed during the simulation. Of 
these, five are BS-BS systems. There appears to be two fairly 
distinct BS binary populations, one consisting of close circu- 
lar orbits and the other with wider eccentric orbits. The BSs 
in circular orbits formed by stable mass-transfer, beginning 
in general when the primary is in the Hertzsprung gap (HG) . 
These still have their primordial companion which evolves 
to become a WD. No instance of a BS formed from Case C 
mass transfer, or as a result of wind accretion, occurs in 
this simulation. For the relatively wide binaries, P > 100 d, 
the distribution appears uniform for e > 0.2. These binaries 
form in hierarchical systems or as a result of an exchange 
interaction. A star is more likely to be exchanged into a wide 
orbit than into an existing hard binary. We also expect that 




3000 3500 4000 

t/Myr 

Figure 9. The evolution of the number of blue stragglers (full 
line), triple systems (dashed line), blue straggler binaries (dash- 
dot line), and RS CVn systems (dotted line), in the M67 simula- 
tion. 
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Figure 10. Distribution of periods and eccentricities for binaries 
formed with a blue straggler component during the M67 simula- 
tion. The open points represent BS-BS binaries, i.e. each compo- 
nent is a BS. 



the eccentricities of newly formed binaries follow a thermal 
distribution (Heggie 1975), which is proportional to the ec- 
centricity. Furthermore, since exchange interactions tend to 
take place in the dense central regions of the cluster it is 
not surprising to find a paucity of wide circular binaries. 
A total of 53 exchange interactions were recorded during 
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L °gl0( T eff/ K ) 

Figure 11. Hcrtzsprung-Russcll diagram for the M67 iV-body simulation at an age of 2 932.5Myr when 3 894 single stars and 3 953 
binaries remain. Main-sequence stars (dots), blue stragglers (open stars), sub-giants, giants and naked helium stars (open circles) and 
white dwarfs (dots) are distinguished. Binary stars are denoted by overlapping symbols appropriate to the stellar type of the components, 
with main-sequence binary components depicted with filled circles and white dwarf binary components as © symbols. The effective 
temperature of a binary is computed according to Hurley & Tout (1998). 



the simulation. Of these, 18 are the result of a direct ex- 
change and the remaining 35 are the result of 12 resonant 
exchanges (Heggie 1975). In a resonant exchange the three 
(or four) stars involved in the interaction form a temporary 
bound system which can then undergo a series of exchanges 
in quick succession before one (or two) of the stars escapes. 

On average half the BSs formed during the simula- 
tion would have done so without the introduction of cluster 
dynamics. We note however that the cluster environment 
was taken into account when choosing the binary parame- 
ters for the starting model, effectively doubling the number 
of BSs expected when compared with population synthe- 
sis run PS6. Consider the 16 BSs present when the simu- 
lation ended, of which three are in binaries. Seven of these 
formed from standard Case A mass transfer in a primordial 
binary: the two components coalesced to leave a single BS. 



One formed as a result of standard Case B mass transfer 
and is in a circular orbit with a 17 d period with its original 
companion, now a WD. Another two formed from Case A 
mass transfer but only after a series of close encounters had 
altered the orbital parameters and caused each system to 
circularize. Without perturbations to their orbits these stars 
would not have been close enough to interact. The remaining 
BSs are the result of collisions within binary systems, three 
after a star was exchanged into a highly eccentric orbit and 
the other three after perturbations to the existing binary 
increased the eccentricity and caused the orbit to become 
chaotic. One of the BS binaries that formed from a dynam- 
ical interaction has a period of 813 d and e = 0.3 and the 
other has P = 6.6 d and e = 0.7. The BS in the wider of these 
has a mass Mi = 1.6 M© and was created after Case A mass 
transfer in a primordial binary lead to coalescence. Later in 
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Figure 12. As Figure [ll] at 3 076.7 Myr. Present in the cluster at this time are 21 blue stragglers, nine of which are in binaries. One 
GB star has a BS binary companion in an eccentric orbit. A CHeB star and a WD in an eccentric binary is also present. There are ten 
cataclysmic variables and ten double-degenerate systems. 



the cluster evolution it underwent an exchange interaction 
with a wide binary to form the existing binary. The com- 
panion is a MS star with M2 = 1.3 Mq and the star that 
left the temporary three-body system was the lightest of the 
three, M3 = 0.7 Mq. The other binary was originally part of 
a three-body hierarchy consisting of an outer component in 
an eccentric orbit about a primordial binary with P — 301 d 
and e = 0.18. Inclination-induced perturbations drove the 
inner eccentricity up to 0.94 at which point one of the inner 
MS components collided with the third MS star, producing 
the BS, which remains bound to its original companion. 



Figures [n], 0, |T|, 0, [l| and |l| show the Hertzsprung- 
Russell diagram of the cluster model at various epochs. The 
spatial distribution at 2 933 Myr and at 4 302 Myr is shown 
in Figures ^ and [Hj] respectively, as the cumulative radial 



profiles in XY- and YZ-planesQ It is evident that the BSs 
are concentrated towards the core of the cluster. 

At T = 4 014 Myr (Figure |l|) there are two super-BSs 
in the cluster. One of these has a mass 2.3 times the then 
turn-off mass of 1.32 Mq. The other has a mass of 2.7Mto 
and is in a binary with another BS. A total of eight super- 
BSs formed during the simulation. Also present in the cluster 
at T — 4 014 Myr are nine cataclysmic variables and three 
double-degenerate systems. There are 19 giants, six of which 
are in binaries. One of the 27 BSs is in a binary with an 
AGB star separated by 4 x 10~ 4 pc. This is about 0.1" at 
the distance of M67, so it would not appear as a BS on a 
cluster CMD (the CCD used by MM J had 0.77"/pixel). 



5 Full colour versions of the Hcrtzsprung-Russell diagrams and 
the spatial distributions, at a larger number of epochs, can be 
found at www.ast.cam.ac.uk/~cat/m67.html. 
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Figure 13. As Figure [TT| at 3 653.4 Myr. Present in the cluster at this time are 29 blue stragglers, the maximum number during the 
simulation, with eight of these in binaries. There is one super-BS, five cataclysmic variables and seven double-degenerate systems. The 
depletion of double-degenerate systems is due to ejection of some close systems and break-up of some wide systems. 



The 2.7M T o super-BS in the model at 4 014 Myr is an 
interesting case. At the beginning of the simulation it is a 
single star with mass Mi = 1.33 Mq . At 3 190 Myr it enters a 
triple system in which it exchanges into the original binary. 
The resulting binary has orbital parameters P = 9 120 d 
and e = 0.52 and survives to 3 740 Myr when it is involved 
in a binary-binary encounter. This leaves the proto-BS in an 
orbit of P = 153 d and e = 0.99 with a companion of mass 
M2 = TOM©. Owing to the large eccentricity the binary 
components collide and merge to form a BS of mass Mi = 
2.33 M Q . At 3 870 Myr this BS is involved in a three-body 
interaction and forms another binary, this time with a star of 
mass Mi = 2.28 Mq which is itself a BS. The BS-BS binary 
has P = 1 350 d and e = 0.67. At 3 990 Myr this binary forms 
a stable four-body hierarchy with yet another binary, the 
original BS collides with one of the stars in the other binary, 
forms the super-BS with Mi = 3.45 Mq, and remains bound 
to the other BS. The fourth body escapes and leaves the 



binary that is observed at 4 014 Myr. Later at 4 080 Myr the 
super-BS and BS collide when the eccentricity of the orbit 
has grown to 0.93 as a result of external perturbations. This 
makes a BS with Mi = 5.7 M Q , i.e. 4.4M T o- The new super- 
BS captures a companion at 4 112 Myr, with which it then 
collides to form a super-BS with Mi = 7.7 Mq ~ 6Mro- 
Soon after it strongly interacts with a hard binary and gets 
ejected from the cluster. 



7 DISCUSSION 

An obvious problem with our M67 model is the lack of BSs 
in wide circular binaries because two such systems, S975 
and S1195, are observed in the real cluster. Leonard (1996) 
suggests these may be the result of mass transfer from a 
core helium-burning (CHeB) primary star. This does not 
seem likely because a star does not get any bigger during 
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Figure 14. As Figure [n] at 3 797.5 Myr. Present in the cluster at this time are 28 blue stragglers with five of these in binaries. One BS 
binary has a GB star companion which is currently filling its Roche-lobe and transferring mass. There is one super-BS, seven cataclysmic 
variables and seven double-degenerate systems. 



CHeB than it did on the GB, so if the primary does not fill 
its Roche-lobe on the GB, i.e. Case B mass transfer, Case C 
mass transfer is very unlikely before the star has reached the 
AGB. If mass transfer does occur when the primary is on the 
AGB then common-envelope evolution can only be avoided 
if the primary has lost enough of its envelope to become the 
less massive star. Wide binaries can form via exchange in 
the cluster but the results show this is unlikely to produce 
circular orbits. So it would seem that a number of wide cir- 
cular BS binaries should be formed by isolated evolutionary 
processes, either from stable Case C mass transfer or wind 
accretion. However, hardly any such binaries are produced 
with the period distributions used in the population syn- 
thesis runs. Perhaps this indicates that a bi-modal period 
distribution is required, or that the peak in the distribution 
should cover a wider range of periods. A more likely solu- 
tion is that the wind-accretion efficiency we have used is too 
low. Our treatment of wind accretion in the binary evolution 



model assumes that the wind velocity Vw is simply Ksc/\/2 
where V CBC is the escape velocity from the surface of the star 
(see Section 2.1 of Hurley, Tout & Pols 2000 with (3 = 0.5). 
The accretion rate is approximately proportional to V„ 4 so 
that a small error in V w has a large effect. We repeated run 
PS6 with the wind velocity reduced by a factor of 2, i.e. 
(3 = 1/8 in Section 2.1 of Hurley, Tout & Pols (2000). This 
is actually more in keeping with observations. The result is 
A BS = 5.1 and N RS = 4.3 at 4.2 Gyr with 23% of the BSs 
in close binaries and 16% in wide binaries (P > lyr). The 
average period of these wide BS binaries is 2 000 d which is 
consistent with the M67 BS binaries. In future A^-body simu- 
lations we will incorporate this variation. Furthermore, if we 
use the weaker criterion of Webbink (1988) to determine the 
onset of dynamical mass transfer at RLOF (see Section 2.6.1 
of Hurley, Tout & Pols 2000) , an additional channel for the 
formation of binary BSs is created. Repeating run PS6 with 
g cr it according to Webbink (1988), together with j3 = 1/8, 



© 2000 RAS, MNRAS 000, 



20 J. R. Hurley, C. A. Tout, S.J. Aarseth and 0. R. Pols 



t = 4013.8 Myr 

N = 1063 N K = 1278 

s b 




3.8 



3.6 



L °gl0( T eff/ K ) 

Figure 15. As Figure [n] at 4 013.8 Myr. Present in the cluster at this time are 27 blue stragglers with six of these in binaries. There 
are two super-BSs and one BS-BS binary. There are nine cataclysmic variables and three double-degenerate systems. 



gives N BS = 5.5 at 4.2 Gyr with 22% of the BSs in close 
binaries and 22% in wide binaries. 

None of the eccentric BS binaries in the model are 
the result of stable mass transfer followed by perturbations 
to the circular orbit. Primarily this is because no BSs are 
formed in wide circular binaries in the first place, but even if 
they were the eccentricity induced would be small and tidal 
forces would quickly return the orbit to circularity. Rasio & 
Heggie (1995) show that the induced eccentricity ef varies 
as 



ef oc 



-5/2 



(24) 



where r p is the separation at periastron. Therefore an eccen- 
tric orbit is more likely to be altered by a close encounter. 
Assuming a stellar number density of 20 pc -3 and a velocity 
dispersion of 0.40 km s" 1 in the core of M67, Leonard (1996) 
uses the result of Rasio & Heggie (1995) to derive a mean in- 



duced eccentricity of e ~ 10 for wide, P = 10 d, circulat- 
or bits. 

Not surprisingly most of the eccentric BS binaries are 
the result of exchange interactions and none are produced 
from tidal capture. The exchange timescale can be expressed 
as 



(25) 



n b E V rc 



where rib is the binary number density, V la \ is the relative 
speed of the third body and the binary centre-of-mass, and 
E is the exchange cross-section. Consider the likelihood of a 
BS of mass Ms = 2.0 Mq being exchanged into a binary with 
Mi = 1.0 M and M 2 = 0.5 M Q . From the results of binary- 
single-star scattering experiments performed by Heggie, Hut 
& McMillan (1996) the exchange cross-section to displace 
Mi is 13 300 a AU 2 and to replace M 2 is 35 400 a AU 2 . Taking 
«b = 20 pc -3 , Vrci = lkms -1 and a — 100 AU gives r cx ~ 
540 Myr which is an order of magnitude less than the age 
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Figure 16. As Figure |ll] at 4 302.1 Myr (the end of the simulation). Present in the cluster at this time are 16 blue stragglers with three 
of these in binaries. There are six cataclysmic variables and three double-degenerate systems. 



of M67. On the other hand, the tidal capture timescale for 
the BS is about 10 4 Myr (Press & Teukolsky 1977). This is 
confirmed by Portegies Zwart et al. (1997) who find that in 
a cluster core with log (n/pc -3 ) = 3.92 the low encounter 
rate means that tidal capture is rare. We must however be 
careful with the definition of tidal capture because some 
BS binaries in the M67 model do form from bound triple 
systems which themselves formed from a binary-single-star 
interaction. 

Even though during the M67 simulation two BSs, 
formed from a collision within a triple system, are found 
in short-period eccentric orbits it is hard to see how the 
cluster dynamics can produce significant numbers of these 
binaries. The exchange cross-section is proportional to the 
binary separation (Heggie, Hut & McMillan 1996) so they 
are unlikely to be formed in this way (see Figure |ro| ) . Al- 
lowing for the chance formation of a tidal capture binary, 
the calculations of Portegies Zwart et al. (1997) show that 



capture binaries in a dense cluster core are generally close 
with a < 10 R0 but that 60% of these circularize during the 
formation process. So for every BS observed in a close ec- 
centric orbit at least one should be found in a close circular 
orbit. BSs in close circular orbits can be formed from Case B 
mass transfer but these are even less likely to have an ec- 
centricity induced than wide circular binaries. A possibility 
that has yet to be considered and does not rely on clus- 
ter dynamics is formation via common-envelope evolution 
when unstable Case B mass transfer occurs. The processes 
involved in common-envelope evolution are very uncertain 
(Iben & Livio 1993) and there is no real reason to assume 
that a binary emerging from this phase should be circular. 
Therefore the binary S1284 observed in M67 could be pro- 
duced in this way or possibly by a collision within a triple 
system. 

The number of blue stragglers produced by the M67 
model is in good agreement with the observations. Figure 
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Figure 17. Cumulative radial profiles in the XY- and YZ-planes for the population shown in Figure |ll| The tidal radius is 16.3 pc 
(vertical dashed line), the half-mass radius is 3.2 pc and the half-mass radius of the blue straggler stars is 0.8 pc. Note that stars are 
not actually removed from the simulation until they are at a distance greater than two tidal radii from the cluster centre. Both blue 
stragglers and giants congregate towards the centre. The four RS CVn systems present are plotted as filled stars. 



shows that the cluster environment is very effective at in- 
creasing the relative number of the BS population. If it had 
been possible to start a full simulation from T = then 
hopefully the relative number of BSs would increase grad- 
ually with time to reach the M67 point, rather than the 
rapid increase produced by the semi-direct simulation (al- 
though the log-scale in Figure [j] distorts the actual range of 
the simulation). As expected the BSs formed by a variety 
of paths, Case A and Case B mass transfer in binaries with 
orbital parameters unaffected by the cluster environment, 
Case A mass transfer after perturbation-induced circular- 
ization, and collisions in highly eccentric binaries. 

In population synthesis without cluster dynamics the 
dominant formation mechanism is Case A mass transfer. 
However, Mathys (1991) has observed that M67 BSs rotate 
slower than average for MS stars. As noted by Leonard & 
Linnell (1992), tidal synchronization in a close binary should 



cause BSs produced by Case A mass transfer to be rapid ro- 
tators. Magnetic braking will not be effective at reducing the 
rotation rate because MS stars with Mq do not have 

convective envelopes. It is uncertain whether BSs resulting 
from collisions between two MS stars will be rapid rotators 
but it is interesting that 50% of the BSs in the TV-body 
model came from direct collisions in eccentric binaries. The 
progenitor stars will not have been affected by standard tidal 
circularization in this case but it is unclear what the effect of 
any angular momentum exchange during chaotic motions of 
the orbit will have on the stellar spins. Our simulation did 
not form any BSs from hyperbolic collisions between two 
single stars, in agreement with the collision timescale pre- 
dicted by Press & Teukolsky (1977), which is about 10 6 Myr 
for a 1 Mq star in the core of M67. 

An encouraging number of wide eccentric BS binaries 
are formed during our simulation. Super-BSs are also made, 
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Figure 18. Cumulative radial profiles in the XY- and YZ-planes for the population shown in Figure The tidal radius is fO.f pc 
(vertical dashed line) and the blue straggler half-mass radius is 0.9 pc. The RS CVn system present is plotted as a filled star. 



as well as some short period binary BSs in eccentric and 
circular orbits. So the model contains formation paths for 
all the BSs observed in M67, except wide circular BS bina- 
ries but this can be rectified as discussed above. Another 
discrepancy is that the frequency of model BSs found in bi- 
naries is too low. The BS half-mass radius in the model is 
smaller than the half-mass radius of the MS stars, in quali- 
tative agreement with the observations, but is itself a factor 
of 2 less than is observed for M67. Also too many potential 
RS CVn systems are broken-up in the cluster core which sug- 
gests that the population within the half-mass radius is too 
centrally condensed. Observations of M67 suggest that the 
number density of stars is about n = 40 pc -3 within a core of 
radius lpc, while the model at 4 300Myr has n = 150 pc~ 3 
within lpc of the cluster centre and n = 34 pc inside the 
half-mass radius. A decrease in the central concentration im- 
plies that wide binaries formed from exchange interactions 
have longer lifetimes. 



The use of a rather unusual tidal field in the M67 model, 
prompted by the observed boundary of 10 pc for the clus- 
ter stars, deserves further discussion. Preferential escape of 
low-mass stars from the model owing to two-body encoun- 
ters shows that the current observed mass of M67, for stars 
with masses greater than 0.5 M©, may be close to the actual 
cluster mass, and therefore a standard tidal field may well 
apply. The case for this is strengthened when we consider 
that the observed boundary of the cluster members is only a 
lower limit for the tidal radius. However, it is interesting that 
v G = 350 kms -1 is not ruled out by the model of the Milky 
Way halo presented by Wilkinson & Evans (1999). Addition- 
ally, Baumgardt (1998) has shown that it is the tidal radius 
at perigalacticon of an eccentric orbit that determines the 
dissolution of a star cluster. The main effect of a stronger 
tidal field is to drive the cluster evolution at a higher rate. 
Therefore it is possible that use of a standard tidal field 
would cause the peak in N B s that we see in Figure to 
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occur at a later time, closer to the age of M67. It is also pos- 
sible that as the core would not be as dynamically evolved it 
would be less dense, leading to a greater population of wide 
binaries and RS CVn systems. We should stress that these 
points are purely conjecture and that the tidal field requires 
close attention in future simulations. In particular, work is 
currently underway to implement a time-varying tidal field 
in NB0DY4 which will enable eccentric cluster orbits to be 
followed (Wilkinson & Hurley 2001). 

Ideally it would be desirable to perform more simula- 
tions and investigate the effects of varying parameters such 
as the central concentration in the starting model. Many fac- 
tors are uncertain, both in the model and the observations. 
For example, the derived escape rate may not be correct for 
larger TV so the size of the starting model, in terms of both 
star number and length scale, may not be relevant to the 
conditions of M67 at birth. In all of this it would be helpful 
if we could begin a complete simulation from initial con- 
ditions but, considering that a single semi-direct simulation 
took a month to perform, the capability of the current hard- 
ware is already pushed to its limit. Future improvements in 
computing efficiency, such as the availability of GRAPE-6 
(Makino 1999), will enable a leap forward in the field of star 
cluster modelling. The ability to perform direct simulations 
of clusters like M67 comfortably will decrease substantially 
the uncertainties involved, and allow observed and model 
parameters to be matched iteratively. 

Finally we also need an explanation for the two sub- 
subgiant binaries observed in M67, S1063 and S1113. These 
lie below the base of the GB (BGB) in the CMD and to 
the right of the MS, much further displaced than the binary 
MS. S1063 has orbital parameters P = 18.4 d and e = 0.217 
while S1113 is circular with P = 2.8 d. HG stars respond to 
changes in mass on a thermal timescale (see Section 2.7 of 
Hurley, Pols & Tout 2000) so that as they lose mass they will 
evolve below the HG of the cluster isochrone. Mass transfer 
is consistent with the orbit of S1113 but not the eccentric 
orbit of S1063, unless the observed eccentricity could be for 
the outer orbit of a triple system. The population synthe- 
sis run PS6 produces 30 binaries in the subsubgiant area at 
4.2 Gyr per 500 000 evolved. These are all circular, have a 
HG primary of average mass 1 Mq, and an average period of 
2 d. So if M67 originally had 15 000 binaries then one binary 
with these orbital parameters can be expected, i.e. S1113. 

Possible evolution paths for S1063 are harder to find. 
The sub-luminous star may have formed in the same way as 
for SI 113 and then been exchanged into an eccentric binary. 
However it is unlikely that two HG primaries are transferring 
mass at the same time and the exchange timescale for a 
short-period binary is long as has already been discussed. 
Another explanation is that a MS star with M < Mro is 
evolving off the MS before it should. Consider a primordial 
binary with Mi = 1.4 M , M 2 = 0.4 M Q , P = 10 d and 
e = 0.9. At T = 2 500 Myr the more massive star begins 
to transfer mass to its companion. The normal MS lifetime 
of the primary is 3 370 Myr but due to the mass transfer it 
ages as it loses mass and at 2 700 Myr, when Mi = 1.2 Mq, 
its effective MS lifetime is 4 200 Myr. If for some reason the 
mass transfer were halted at this point then the star would 
evolve off the MS at 4 200 Myr into the subsubgiant region. 
Exchange of the star into a wider binary is one way to halt 
the mass transfer and to explain the observed parameters of 



S1063. In run PS6 there are 600 likely exchange targets per 
500 000 binaries but the window for exchange is small, as it 
must occur when Tms — 4 200 Myr, and so is the exchange 
cross-section. The exact nature of S1063 remains a mystery. 



8 CONCLUSIONS 

In order to complement increases in computing speed we 
have made a substantial effort to improve the treatment of 
physical processes within the cluster environment. Our TV- 
body code includes detailed modelling of stellar and binary 
evolution, thus allowing us to test directly the influence of 
the interaction between stellar evolution and gravitational 
encounters on the evolution of a star cluster. Some simi- 
lar advances have also been made in the work of Portegies 
Zwart et al (1998, 2000) although their models still only in- 
clude Population I evolution and are therefore not suited to 
studying globular clusters. In particular, certain aspects of 
binary evolution, such as tidal circularization, are not mod- 
elled in as much detail as this work. 

NB0DY4 is extremely useful for simulating cluster pop- 
ulations. As a first application we have modelled the blue 
straggler population in M67. We could just as easily have 
looked at many other aspects of cluster evolution, such as 
the production of CVs or development of mass segregation, 
but such topics will be the subject of future work. In par- 
ticular, the availability of GRAPE-6 will enable us to model 
small globular clusters directly. 

We have shown that binary evolution alone cannot ac- 
count for the numbers of observed blue stragglers in open 
clusters, or the binary properties of these blue stragglers, 
when a realistic separation distribution is assumed. The in- 
fluence of the cluster environment can effectively double the 
number of blue stragglers produced, leading to good agree- 
ment with the observations. Our TV-body model of M67 
demonstrates that blue stragglers are most likely generated 
by a variety of processes and in particular we find formation 
paths for all the BSs observed in M67. We also find that, 
among the possibilities we consider, the primordial binary 
population is best represented by a log-normal distribution 
of separations peaked at 10 AU and binary masses chosen 
from the mass function of Kroupa, Tout & Gilmore (1991) 
in combination with a uniform mass-ratio distribution. 

We quantify the escape rate of stars from a clus- 
ter subject to the tidal field of our Galaxy as M = 
— 0.3M / (i r h l°gi()0-4TV) . This enables us to provide a method 
by which the initial cluster mass and radius corresponding 
to current values can be determined. So we were able to 
model M67 by a semi-direct method even with current com- 
putational limitations. 
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